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Abstract 

The invariant formulation of CP violation involves the generation of sets of invari- 
ant constraints for CP conservation, the manipulation of their expressions and the 
identification of complete and minimal subsets of such constraints. In this paper we 
' present a collection of subroutines to deal with these three tasks in a fast, reliable and 

\Q . systematic way, with examples for the leptonic sector. 
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2h ■ 1 Introduction 



The observed CP violation in the K° — K° system [|J is related to the presence of complex 
phases in the mixing matrix describing the quark gauge couplings in the mass-eigenstate 
basis in the Standard Model (SM) Q. However, not all the phases in the mixing matrices 
are CP violating. Some of them can be eliminated redefining the fermion phases. In the 
SM with three quark generations, the 6 phases of the 3x3 unitary mixing matrix V 
reduce to 1 after an appropriate field redefinition. For this case it has been identified a 
quantity invariant under weak quark basis transformations whose vanishing characterizes 
CP conservation ||: 

det [M u Ml,M d M ] d ] = -2i(m 2 - m 2 c )(m 2 - ml)(m 2 c - m£)(mf - m 2 s ) 

x (ml-m 2 d )(m 2 s -m 2 d )Im(V ud V c * d V cs V: s ) = 0, (1) 

where m; is the mass of the quark i and Vij is the ij entry of the Cabibbo-Kobayashi- 
Maskawa (CKM) matrix H. This invariant formulation of CP violation makes apparent 
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the necessary and sufficient conditions for CP conservation, giving the size of CP violation. 
It also allows to decide if CP is conserved in any weak basis, motivating model building 
and helping to understand eventually the origin of CP violation if a definite model (weak 
basis) if physically distinguished. 

In general, i.e. in extended models with extra generations and/or vector-like quarks 
and/or right-handed (Majorana) neutrinos, it is difficult to find such a minimal set of 
necessary and sufficient invariant conditions. In this paper we present a collection of 
subroutines for Mathematica pj] which carries out this task in a fast and systematic way. 
They are used (i) to calculate products of mass (sub)matrices giving rise to expressions 
invariant under weak fermion basis transformations, (it) to do the symbolic manipulations 
and to solve the constraints resulting from equating these invariants to zero, and (iii) to 
verify that the set of necessary conditions for CP conservation generated in this way is 
sufficient and minimal (this is to say that this set does not contain any smaller subset of 
sufficient conditions). Although these subroutines are written for the leptons including 
Majorana neutrinos, they can be easily adapted to analyze CP violation in the quark 
sector H [|. They can be also generalized to study extended gauge theories as left- 
right models, and/or (supersymmetric) scalar gauge couplings. On the other hand, the 
computer calculations use to be more reliable than the manual ones. For example, we 
have revised the two more complicated cases studied in Ref. J5J. In both cases we find CP 
violating solutions which were overlooked there. In Section 2 we comment on the main 
theoretical features and in Section 3 we describe the subroutines and apply them to these 
two cases. 



2 Theory and general features 

Let us assume the SM with hl left-handed and ur right-handed neutrinos and define 

where v c = i^ 2 v* and ipL,R are til + n R component vectors describing the til left-handed 
and riR right-handed neutrinos. Then the mass term reads 

- C mass = IlMJr + -tPlM^r + h.c, (3) 

where Il,r are component vectors describing the charged leptons and Mi an arbi- 
trary til x til complex matrix, whereas 

riL n R 

_ n L { (Ml 5&\ 
Mv ~ n R { \ mIMr) W 
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is an arbitrary (til + rift) x (til + ur) complex symmetric matrix. This Lagrangian 
is invariant under a CP transformation leaving the gauge SU(2)l x U(1)y interactions 
unchanged 

h - L/lC£ , Z fl - C^C7£ , ^ - £/lCV£ , - C^C^ , (5) 

where C is the Dirac charge conjugation matrix and Ul, U l R (U R ) are ul x til (ur x hr) 
unitary matrices, if 

U^Ur = Ml , u[M L Ui = M* L , uIMdUr- = M* D , U V R T M R U» R = M* R . (6) 

Eqs. (||) are necessary and sufficient conditions for CP conservation. From these equations, 
one can write necessary conditions for CP conservation which do not require the knowledge 
of the unitary matrices involved in the definition of the CP transformation. Thus, the 
products of mass matrices Mi, Ml, Mo, Mr are classified in three classes G\ = Gl, 
G2 = G l R , G3 = G R , depending under which unitary matrix (Ul, U r , U r , respectively) 
transform, 

Gi = {A L i = M l M},A L2 = M L M\, A L3 ee M D M D , A Li A Lj , M L M?M?Ml 

m l m* d mIm ] l , m l m* d m r m ] d , m d m r mIm\, M d M r M r M ] d , . . .}, 
g 2 = {a 1 = m}m 1 ,aIm}m l mIm 1 ,mJm d mI ) m 1 ,...}, 

G 3 = {A ul ee M^Md, A u2 ee M R M R , A vi A„j , M^MiM} M d , M ] d M l m[m d , 

m ] d m l m* d m r , m r mIm\m d ,m r mIm* d Mr, . . .}. (7) 

The trace and determinant of any of such products or their sums within the same class do 
not change under unitary transformations due to their cyclic and factorization properties, 
respectively. Then CP conservation (Eq. (||)) implies that the trace and the determinant 
of any linear combination g of elements of the same class is real, 

Im tr g = , Im det g = 0. (8) 

To find which subsets of these conditions are also sufficient has to be worked out case by 
case. In practice, we are interested in identifying minimal subsets with these properties. 
We proceed as follows. 

2.1 Generation of invariant products with the same transformation prop- 
erties 

The elements of Gi in Eq. ([7]) must be generated using the mass matrices 

{M h Mj, Ml, M\, M L , Ml, M D ,M T D , M* D ,M^ D , Mr, M r }. (9) 
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(Note that Ml <R are symmetric.) A sequence S of order n is an allowed product of n of 
such matrices. A product of two sequences S a S\, is allowed if S' a = V a S a V§ , S' b = V&»S&V^ 
and V' a = Vb- Hence the sequences of order 2 are 

{M/M/, Mf M* , Mf M^, M?M^,M?M?, M}M h M}M L , m}m d , 
M L M?,M L Ml M L M* D ,MlM h m\M u mIm d , M D Mjy, M D M R , 

mIm^mImI, m£m* d , m* d mI, m* d m r , m ] d m u mIm l , M ] d M d , 
m r mI,m r m r ,m r m£,m r m r }. (10) 

A sequence S, with S' = VSV'^ is an element of £71,2,3 in Eq. ([?]) if V = V = 
Ul, U l R , U R , respectively. Then Eq. @ implies that the order of the elements of G{ is even 
and that the sequences in Gi can be constructed with the biproducts in Eq. (|l0|). These 
elements are generated order by order. The elements of order n are the allowed products 
of the sequences of order n — 2 which begin with Mi, Ml or Mb for G±, Mj for G2 and 
Mp or M R for G3 (see Eq. @), times the biproducts in Eq. (|i"o|) which end with Mj, 
m\ or M ] D for G x , M t for G 2 and M D or M R for G 3 . 

2.2 Solution of invariant constraints 

Once the elements of Gi are generated, we have to solve the nontrivial conditions in 
Eq. (||) up to a given order. This means to find the relations among the parameters 
fixing Mi^l,d.r implied by these conditions. The order is increased until there is no CP 
violating solution. In the examples discussed later we need to consider conditions involving 
sequences of order 4, 6 and 8. These constraints are easier to solve if the mass matrices 
are conveniently parametrized. Under a change of weak basis, 

M{ = W L MtW 1 ^ , M' L = W L M L W[, 

M' D = W L M D Wtf , M' R = WgM R Wg . (11) 

Hence, as the constraints do not depend on the choice of weak basis, we can choose the 
unitary matrices Wl, W r u appropriately and assume Mi and M R diagonal with nonnega- 
tive real elements, whereas Ml is complex and symmetric and Mjj complex and arbitrary. 
The number and difficulty of the equations grow with the order of the sequences. Gener- 
ically Mathematica uses too much memory and time to solve simultaneously more than 
5 conditions. In this case one solves the equations of lowest order first and inserts the 
solutions in the remaining equations, which can then be solved more easily. 
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2.3 Minimal subset of invariant constraints 



The last subroutine verifies if a solution of a set of conditions in Eq. (|8|) is CP conserving. A 
solution s is a set of relations among the parameters fixing the mass matrices Mil d R- Let 

(s) 

M[£ D R be the mass matrices satisfying these relations. We generate random numbers for 

' ' ' ~ (s) 

their independent parameters denoting the corresponding mass matrices M^ LD R . These 
matrices conserve CP if it exists a unitary matrix preserving M, and such that 



(12) 



If this matrix can be constructed, we say that s conserves CP. It is highly improbable 

~ (s) 

that there is a relation among the random parameters in M^ £ D R such that the matrices 

(s) ~ (s) 

M^£ D R violate CP but £ D R do not. On the other hand, if Eq. ( jl2| ) is not satisfied, 
the solution s is said to violate CP. 

Let us drop the superscripts from now on. The existence of U is determined in two 
steps. Step 1 guarantees that if U exists, it is diagonal, and step 2 checks if there is 



a diagonal unitary matrix satisfying Eq. (12). In the convenient basis and for a given 
solution s, U is block diagonal with diagonal submatrices of dimension the multiplicity rrii 
and /ij of the eigenvalues of M; and Mr respectively, 



U 



{m\ x mi) 






(ni2 x m 2 ) 



\ 



(/XI X fli) 






(/U 2 X fj, 2 ) 



V 



(13) 



T2|) is fulfilled if 
(14) 



with J2 m i = n L , E Mi = n R- 

If there is no degeneracy U is diagonal and we go to step 2: Eq. 
Vi > j, 

= arg Ma - arg My + arg Mjj - arg M^ = 0, tt. 

If any M y - = 0, A y - is zero for argMjj (= arg Mjj) is conventionally defined by Eqs. 
(|l4|). On the other hand, if some diagonal elements = • • • = Mi k i k = 0, we consider 

argMj^, . . . , argMj fc j fe arbitrary but fixed. We regard Eqs. (|l4|) as a system of equations 
in the variables Si = argMjj. Then CP is conserved if and only if the full system of Eqs. 
(Ilil) is compatible. 

If there is degeneracy we can fix the weak basis to ensure that U is diagonal and 
then go to step 2. This involves the most delicate casuistry, and it has to be done and 
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programmed for each dimension tir. For small dimensions it can be easier to do it 
by hand. Let us discuss the simplest case: ul = 1, tlr = 2, and the 2x2 matrix Mr 
diagonal, real, nonnegative and degenerate. Then U has the general form 



"={ o ul)- < 15 > 

If Mr is identically zero, Ur is an arbitrary unitary matrix. There are different ways 
to fix the weak basis and guarantee that Ur is diagonal. We choose to diagonalize 
M D M D . 

— If M D Mp is nondegenerate, we transform M accordingly and go to step 2. 

— If MqMd is degenerate, it is identically zero because Mo is a 1 x 2 matrix. 
Then Mb is also identically zero and there is an extra flavour symmetry, and 
Eq. ( |T2] ) is satisfied taking Ur equal to the identity. 

If Mr is nonzero, Ur is a real orthogonal matrix. In this case we choose to diagonalize 
(Im Md) t (Im Md) or (Re Md) t (Re Mp) if the former is degenerate. 

— If one of them is nondegenerate, we transform M accordingly and go to step 
2. (For convenience only, if Ml = we also require MpM]^ to be real and 
positive. This fixes 5l- If MdMJ) = there is a flavour symmetry.) 

— If both are degenerate, as above they are identically zero. Mo is also zero and 



there is again an extra flavour symmetry. Similarly Eq. (12) is fulfilled with 
Ur equal to the identity. 



3 Examples 

In this section we present a description of the different subroutines in the package CPlep. 
Some of them require the subroutines DiagonalizeH and DiagonalizeS defined in the 
package Diagon |J. DiagonalizeH is a subroutine which diagonalizes hermitian ma- 
trices, returning orthonormal eigenvectors even in the case of degenerate eigenvalues^. 
DiagonalizeS diagonalizes general complex symmetric matrices M s with a congruent 
transformation M s — > UM S U T = D. Let us show how CPlep works in two examples. 
(Definitions are given in the Appendix.) 

1 The built-in function Eigensystem may give non-orthogonal eigenvectors 
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3.1 Case til — 1 5 tir — 2 
The list 



ll={a[10] ,a[ll] ,a[12] ,a[13] ,a[20] ,a[23] ,a[30] ,a[31] ,a[32] ,a[33] , 
a [40] ,a[43]> 

is the computer representation of the mass matrices in Eq. (||). We define an adequate 
output format for a[i] so that the expressions can be more easily understood. The rules 
for the generation of invariants are embodied in the definition of an associative function 
f , with an arbitrary number of arguments a[ii] ,. . . ,a[z' m ] ordered as the corresponding 
matrix product a[«i] • • -a[i m ] . f is taken to be zero if any of the products a[ii] a [22] , 
... ,a[i m _i] a[i m ] is not allowed (see Eq. (|l0|)). Then the function inv checks if the 
product f [a[ii] ,. . . ,a[i m ]] transforms with a matrix and its adjoint, and returns again 
its argument if it does or if it does not. The lists obtained in this way are multiplied 
with Outer. 

«cplep.m; 

12 [1] =Union [Flatten [Outer [f , cllini , 11] ] ] ; 
12 [2] =Union [Flatten [Outer [f , cl2ini , 11] ] ] ; 

12 [3] =Union [Flatten [Outer [f , cl3ini , 11] ] ] ; 
13 [1] =Union [Flatten [Outer [f , 12 [1] , 11] ] ] ; 

13 [2] =Union [Flatten [Outer [f , 12 [2] , 11] ] ] ; 

16 [3] =Union [Flatten [Outer [f , 15 [3] , 11] ] ] ; 

14 [1] =Union [inv [14 [1] ] ] ; 
i4 [2] =Union [inv [14 [2] ] ] ; 
i4 [3] =Union [inv [14 [3] ] ] ; 
i6 [1] =Union [inv [i6 [1] ] ] ; 
i6 [2] =Union [inv [i6 [2] ] ] ; 
i6 [3] =Union [inv [i6 [3] ] ] ; 

In this example it is necessary to consider only sequences up to order 6. In order to obtain 
the invariant constraints for til = 1, nji = 2, we introduce the convenient parametrization 
for Ml, ML, MD, MR and define Mn, the neutrino mass matrix in Eq. (Q). Random numerical 
values for the parameters are also generated for later use. 

Ml:={{el}};ML={{nl}};MD:={-[Re[n2]+I i m [ n 2] , Re [n3] +1 Im[n3]»; 
MR:={ {n4,0},{0,n5} }; Mn=Transpose [Join [Transpose [Join [ML, 
Transpose [MD] ] ] , Transpose [Join [MD , MR] ] ] ] ; 
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evl={el->Random [] ,nl->Random [] ,Re [n2] ->Random [Real , {-1 , 1}] , 
Im[n2] ->Random [Real, {-1,1}] ,Re [n3]->Random [Real ,{-1,1}] , 
Im [n3] ->Random [Real , {-1 , 1}] ,n4->Random [] ,n5->Random[] } 

Now to identify those constraints in Eq. (||) which are not identically zero, we define the 
functions Looktrace and Lookdet. Their first argument 1 is a list of sequences, whereas 
the other arguments are the mass matrices with random entries. These are required to 
check if the constraint is identically zero. Looktrace looks for the elements g of 1 which 
satisfy Im tr g ^ 0. Lookdet first finds all non-hermitian combinations g — g' and then, 
returns those with Im det(g-g') ^ 0. Finally Newecs calculates the non-trivial constraints 
returned by Looktrace and Lookdet as a function of the convenient parameters (see the 
Appendix for a detailed description of its arguments). Thus 

Looktrace [i4[l] ,M1/ . evl ,ML/ . evl ,MD/ . evl ,MR/ . evl] ; 
Newecs [{},%, Ml, ML, MD, MR, {el, nl,n2,n3,n4,n5},{n2,n3}] 
ecs={ G / [[3,2]]}; 

returns 

Imtr (MdM^M^mI) = 2rai(n 4 Imn 2 Ren 2 + n 5 Imn 3 Rera 3 ). (16) 

The traces of the other elements of i4 [1] , i4 [2] and i4 [3] are proportional to the trace 
in Eq. (|l6|). Similarly, 

Looktrace [i6 [1] ,M1/ . evl ,ML/ . evl ,MD/ . evl ,MR/ . evl] ; 
Newecs [ecs , °/„ , Ml , ML , MD , MR , {el , nl , n2 , n3 , n4 , n5} , {n2 , n3}] 
AppendTo [ecs , °/„ [ [3 , 2] ] ] ; 

gives 

Imtr (MdM^MrM^MdMI) = 2n 1 (nl Imn 2 Rera 2 + nl Imn 3 Rera 3 ). (17) 

The traces of the remaining sequences in i6 [1] , i6 [2] and i6 [3] do not give new equa- 
tions. Finally 

Lookdet [i4 [3] , Ml/ . evl , ML/ . evl , MD/ . evl , MR/ . evl] ; 

Newecs [ecs , % , Ml , ML , MD , MR , {el , nl , n2 , n3 , n4 , n5} , {n2 , n3} , det->True] 

AppendTo [ecs , °/ [ [3 , 3] ] ] ; 

returns 

Im det{M ] D M D M R M R - M R M^M^M R ) = 2ra 4 n 5 (n 4 - n 5 ) (n 4 + n 5 ) 

x (Imra 2 Ren 3 — Imn 3 Rere 2 ) 

x (Imn 2 Imn 3 + Rere 2 Ren 3 ). (18) 
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When we consider that the set of constraints may be complete, we use Reduce to find 



all the solutions. In this case we try with Eqs. (06-18) equal to zero. Reduce gives a 
lot of redundant solutions, often repeated, for instance nl==0 && n4==0 and nl==0 && 
n4==0 && n5==0, and some inconsistent solutions, like Re[n3]==I Re [n2] . The function 
Eliminatesols gets rid of the redundant solutions; whereas the inconsistent solutions are 
eliminated by inspection. 

solO=Reduce [ecs=={0 ,0,0}]; 
«Algebra/ReIm . m ; 

el/: Positive [el] = True;nl/: Positive[nl] = True; 
n4/ : Positive [n4] = True;n5/: Positive [n5] = True; 
soll=Eliminatesols [solO/ . _ ! =0->True] ; 

sol2=soll/ . { (Re [_] ==Complex [_ , _] x_/ ; Im [x] ==0) ->False , 
(Im [_] ==Complex [_ , _] x_/ ; Im[x] ==0) ->False} ; 

In this case, Reduce gives us 198 solutions (solO). Many of them are redundant and 
we need to keep only 32 (soil). 6 of them are inconsistent and are eliminated in sol2. 
To verify whether the former set is complete, we use the subroutines Redl2 and LookCP. 
Redl2 writes the mass matrix Mn in the basis in which U is diagonal, if it exists. LookCP 
then returns True or False depending on whether U exists. 

Although the whole process has been described in the previous section, we want to 
point out a little trick used in LookCP. When trying to solve the system of equations 
Aij = 0,7r in Eqs. (|i~4|), one could consider instead to solve the system of equations 
sinAjj = 0. However, Mathematica does not give correct results in this case. For this 
reason we look for a subset of linearly independent equations in Eqs. (|l~4|), and solve 
these equations equated to zero (it can be done always if the equations are independent, 
because we can always redefine the variables Si = arg Ma conveniently) . Then the solution 
is substituted in the complete set of equations, checking whether they are equal to 0, it. 

We check if any of the 26 solutions in sol2 violates CP using a loop. It is initialized 
with 

m [number_] : =Mn/ . ToRules [sol2 [ [number] ] ] / . evl ; 
i=0; 

and for each solution we run 
++i ; LookCP [Redl2 [m [i] ] ] 

The result is that none of the solutions in sol2 violates CP. So the set of invariant con- 
straints is complete. We could try to take away any of the equations in ecs and repeat 
the same process, but we would find CP violating solutions. So the set is also minimal. 
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3.2 Case ul = 2, ur = 1 

This case is more difficult to solve because more invariant constraints are necessary. As 
we pointed out in Section 2, the computer memory required to solve all the equations 
simultaneously is too big, so we have to use a different approach. First we solve the 
condition 

Im det ( Ml m£ M; Mj - M L MfMfM{) = 0. (19) 

Then for each solution of Eq. ( |I9| ) we proceed as in the former example but restricting 
the general form of the mass matrices to fulfill this particular solution. For this case we 
have found that 



Im tr (M ] D M L M* D M R ) 


= o, 


Im tr (M R M%MlM L MlM D ) 


= o, 


Im tr {MiM} M l M? M[ M^M R Mp) 


= o, 


Im det {Ml M\ Mi M] - M l M?M?mI) 


= o, 


Im ^(MdM^MlMI - M l M* d MdMI) 


= o, 


Im det (Md M^ d Mi M] - M l M* d mIm ] l ) 


= o, 


Im det(M D M j D MiMj - MiM\m l m\) 


= o, 


Im det(M D MlMiMj - M l M^MJm\) 


= o, 


Im det{M L M* D M R M ] D - MiM}MiM}) 


= 



form a complete set of constraints. This set is also minimal by construction, as long as 
Eq. (|19|) is included. 
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A Appendix 



(* 

This package is available by anonymous ftp at deneb.ugr.es 

in directory pub/packages 

*) 

BeginPackage ["CPlep' " , "Diagon' "] 
(* 

This package needs the functions DiagonalizeS , DiagonalizeH 
defined in the package Diagon, available by anonymous ftp at 
deneb.ugr.es 

DiagonalizeS [m] gives a list {v,u} where m is a square matrix, 
v is a row vector and u is a matrix fulfilling 
u . m . Transpose [u] =DiagonalMatr ix [v] 

DiagonalizeH [m] is a modification of Eigensystem which gives a 
list {v,u}, with the eigenvalues and eigenvectors of the 
hermitian matrix m, such that 

u. m. Transpose [Conjugate [u] ] =DiagonalMatrix [v] . The list of 
eigenvalues v is ordered by increasing absolute value 
*) 

cllini : : usage = 

"cllini is a list 

the class Gl" 
cl2ini : : usage = 

"cl2ini is a list 

the class G2" 
cl3ini : : usage = 

"cl3ini is a list 

the class G3" 
cllfin: : usage = 

"cllfin is a list 

the class Gl" 
cl2fin: : usage = 

"cl2fin is a list 

11 



with the possible first elements of a product in 



with the possible first elements of a product in 



with the possible first elements of a product in 



with the possible last elements of a product in 



with the possible last elements of a product in 



the class G2" 
cl3fin: : usage = 

"cl3fin is a list with the possible last elements of a product in 
the class G3" 
11 : : usage = 

"11 is a list with all the possible elements of a product" 
f : : usage = 

"f [expr] is if expr is not an allowed product" 
inv: : usage = 

"inv[expr] returns expr if it is an invariant product, or if not" 

Looktrace :: usage = 
"Looktrace [If ,M1,ML,MD,MR] returns the elements in a list If of 
f-products with trace nonzero when evaluated with the numerical 
matrices Ml, ML, MD, MR. The output is in the form of Dot products" 

Lookdet : : usage = 
"Lookdet [If ,M1,ML,MD,MR] calculates all possible non-hermitian 
combinations x-y of elements of a list If of f-products and then 
returns those with nonzero determinant when evaluated with the 
numerical matrices Ml, ML, MD, MR. The output is in the form of Dot 
products" 

Newecs : : usage = 

"Newecs [oldecs , Id, Ml ,ML,MD,MR, vars , comp] looks for equations not 
present in oldecs, generated with a list Id of Dot products and the 
matrices Ml, ML, MD, MR with variables vars, of which comp are complex. 
It returns a list containing four lists: the list oldecs, possibly 
improved with simpler equations, the list of elements of Id that 
improve oldecs, a list of new equations not present in oldecs and a 
list of elements of Id that generate the new equations" 
det : : usage = 

"det is an option for Newecs with default value False which select 

to calculate the trace or the determinant" 
Eliminatesols :: usage = 

"Eliminatesols [sols] eliminates redundant solutions in sols" 
Redl2 : : usage = 

"Redl2[M] reduces a neutrino mass matrix M in the case nL=l,nR=2 by 
changes of basis to a form in which CP is a diagonal matrix if it 
exists" 
Red21 : : usage = 

"Red21[M,Me] reduces a neutrino mass matrix M in the case nL=2,nR=l 
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by changes of basis to a form in which CP is a diagonal matrix if 
exists. Me is the charged lepton mass matrix" 
LookCP: : usage = 

"LookCP[M] tests if a diagonal unitary matrix U exists such that 
U. M. Transpose [U] =Con jugate [M] and yields True or False" 

Begin [" 'Private' "] 

(* Definitions for internal use only *) 

(* No comments needed *) 

trace [x_?MatrixQ] : =Sum [x [ [i , i] ] , 

{i .Length [x] }] / ; Length [x] ==Length [Transpose [x] ] ; 

Adj [m_] : =Transpose [Conjugate [m] ] 

(* Enlarge the dimension of a matrix with the identity *) 

IncDiml [m_] : =Transpose [Prepend [Transpose [Prepend [m , Table [0 , 
{Length [m] }] ] ] , Join [{1} .Table [0 , {Length [m] >] ] ] ] 

IncDim [m_ ,n_] : =Nest [IncDiml ,m,n] 

(* adjoint of a matrix *) 

adl[a[x_]] :=a[10 Floor [x/10] +3-Mod [x, 10] ] ; 
(* Adjoint of an f -product *) 

ad [f [x__] ] : =Apply [f , Reverse [Thread [adl [Apply [f f , f [x] ] ] , f f ] ] ] ; 
ad[0] =0; 

(* Take hermitian elements of a list *) 
auto [l_List] : =Block [{It} , ( 

For [lt={> ; n=l , n<=Length [1] , n++ , If [1 [ [n] ] ==ad [1 [ [n] ] ] , 
AppendTo[lt,l[[n]]]]] ;lt)] ; 
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(* Take non-hermit ian elements of a list *) 
autono [l_List] : =Block [{It} , ( 

For[lt={};n=l,n<=Length[l] ,n++, If [1 [ [n] ] ==ad[l [ [n] ] ] , , 
AppendTo [It , 1 [ [n] ] ] , AppendTo [It , 1 [ [n] ] ] ] ] ; It ) ] ; 

(* Imtr and Imdet are used by Newecs and just substitute the 
symbolic expressions of the mass matrices and calculate the 
imaginary part of the trace/determinant *) 

Imtr [expr_ , x_ , y_ , u_ , v_ , comp_] : =Block [{sust} , ( 

sust={a [10] ->x , a [20] ->y , a [30] ->u , a [40] ->v , 

a[13]->Adj [x] ,a[23]->Adj [y] ,a[33]->Adj [u] ,a[43]->Adj [v] , 

a [11] ->Transpose [x] ,a[31] ->Transpose [u] , 

a [12] ->Conjugate [x] , a[32] ->Conjugate [u] } ; 

Factor [ComplexExpand [Im [trace [expr/ . sust] ], comp] ] )] 

Imdet [expr_ , x_ , y_ , u_ , v_ , comp_] : =Block [{sust} , ( 

sust={a[10] ->x,a[20] ->y,a[30] ->u,a[40] ->v, 

a[13]->Adj [x] ,a[23]->Adj [y] ,a[33]->Adj [u] ,a[43]->Adj [v] , 

a [11] ->Transpose [x] ,a[31] ->Transpose [u] , 

a [12] ->Conjugate [x] , a [32] ->Conjugate [u] } ; 

Factor [ComplexExpand [Im [Det [expr/ . sust] ], comp] ] )] 

(* minors calculates the set of invariant phases \Delta_{ij} of a 
matrix. Used by LookCP *) 

minor [m_ , i_ , j _] : =m [ [i , i] ] +m [ [ j , j ] ] -m [ [i , j ] ] -m [ [ j , i] ] 

minors [m_] : =Flatten [Table [If [m [ [i , j] ] ==0, , minor [m, i , j] , 
minor [m , i , j ] ] , {i , 1 , Length [m] } , { j , i+1 , Length [m] }] ] 

(* vec and linealindep check if a susbset of invariant phases is 
linearly independent. Used by LookCP *) 

vec [expr_ ,d_] : =Table [Coefficient [expr .delta [i] ] , {i ,d}] 

linealindep [exprs_ , d_] : =MemberQ [Complement [Flatten [Minors [Thread [ 
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vec [exprs , d] ] , Length [exprs] ] ] , {0}] , _] 



(* Output format with standard notation *) 



Format [f [x_ , 


y__]]:={x,y}; 








Format [a [10] 


] =ColumnForm [{" 


ii 

> 




1"}, Left, Center] ; 


Format [a [13] 


] =ColumnForm [{" 


+ " 


,"M", 


" 1"}, Left, Center] 


Format [a [12] 


] =ColumnForm [{" 


*»■ 


,"M", 


" 1"}, Left, Center] 


Format [a [11] 


] =ColumnForm [{" 


J II 


,"M", 


" 1"}, Left, Center] 


Format [a [20] 


] =ColumnForm [{" 


11 

> 


"j^n ii 


L"}, Left, Center] ; 


Format [a [23] 


] =ColumnForm [{" 


+ " 


,"M", 


" L"}, Left, Center] 


Format [a [30] 


] =ColumnForm [{" 


II 

> 


hj^ii ii 


D"}, Left, Center] ; 


Format [a [33] 


] =ColumnForm [{" 


+ " 


,"M", 


" D"}, Left, Center] 


Format [a [32] 


] =ColumnForm [{" 


*" 


,"M" , 


" D"}, Left, Center] 


Format [a [31] 


] =ColumnForm [{" 


J 11 


,"M", 


" D"}, Left, Center] 


Format [a [40] 


] =ColumnForm [{" 


11 

> 


nj^ii ii 


R"}, Left, Center] ; 


Format [a [43] 


] =ColumnForm [{" 


+ " 


,"M", 


" R"}, Left, Center] 



(* Definition of f: not allowed producs are *) 



f [a [10] ,a[x. 


J] 


=0/ 


x! = 


= 13; 








f [a [13] ,a[x. 


J] 


=0/ 


(x 


= 10 


fe& X 


=20 && 


x!=30) ; 


f [a [12] ,a[x. 


J] 


=0/ 


x! = 


= 11; 








f [a[ll] ,a[x. 


J] 


=0/ 


(x 


= 12 


fe& X 


=23 && 


x!=32) ; 


f [a [20] ,a[x. 


J] 


=0/ 


(x 


= 12 


fe& X 


=23 && 


x!=32) ; 


f [a [23] ,a[x. 


J] 


=0/ 


(x 


= 10 


fe& X 


=20 && 


x!=30) ; 


f [a [30] ,a[x. 


J] 


=0/ 


(x 


=33 


&& X 


=43); 




f [a [33] ,a[x. 


J] 


=0/ 


(x 


= 10 


&& X 


=20 && 


x!=30) ; 


f [a [32] ,a[x. 


J] 


=0/ 


(x 


=31 


&& X 


= 40) ; 




f [a [31] ,a[x. 


J] 


=0/ 


(x 


= 12 


&& X 


=23 && 


x!=32) ; 


f [a [40] ,a[x. 


J] 


=0/ 


(x 


=33 


&& X 


=43); 




f [a [43] ,a[x. 


J] 


=0/ 


(x 


=31 


&& X 


=40); 




f [f [x__,y_] 


z_] 


:=If [f [y,z] 


==0, 


D,f [x,y 


,z] ,f [x, 



f [x_,f [y_,z__]] :=If [f [x,y]==0,0,f [x,y,z] ,f [x,y,z]] ; 

f [f [x ,y_] ,f [u_,v ]] :=If [f [y,u]==0,0,f [x,y,u,v] ,f [x,y,u,v]] ; 

f [0,x__] :=0;f [x__,0] :=0; 
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(* Definition of inv *) 



inv [f [x_ , y , z_] ] : =If [f [z , x] ==0 , , f [x , y , z] , f [x , y , z] ] ; 

inv[0] :=0; 

Set Attributes [inv, Listable] ; 
(* Definitions of lists *) 

ll={a[10] ,a[13] ,a[12] ,a[ll] ,a[20] ,a[23] ,a[30] ,a[33] ,a[32] ,a[31] , 
a [40] ,a[43]}; 

cllini={a[10] ,a[20] ,a[30]>; 
cl2ini={a[13]}; 
cl3ini={a[33] ,a[43]}; 
cllfin={a[13] ,a[23] ,a[33]>; 
cl2fin={a[10]}; 
cl3fin={a[30] ,a[40]}; 

(* Definitions of Looktrace and Lookdet *) 

Looktrace [l_List ,x_ ,y_ ,u_ , v_] : =Block[{n, sust , lev, It ,1s} , ( 
ls=Apply[Dot,l,l] ; 

sust={a [10] ->x , a [20] ->y , a [30] ->u , a [40] ->v , a [13] ->Adj [x] , 
a[23]->Adj [y] ,a[33]->Adj [u] ,a[43]->Adj [v] , a [11] ->Transpose [x] , 
a [31] ->Transpose [u] ,a[12] ->Conjugate [x] , a [32] ->Conjugate [u] } ; 
lev=ls/ . sust ; 

For [lt=-Q ; n=l , n<=Length [Is] , n++ , 

If [Chop [trace [lev [ [n] ] ] ] ! =Conjugate [Chop [trace [lev [ [n] ] ] ] ] , 

AppendTo[lt,ls[[n]]] ] ]; 

It)]; 

Lookdet [l_List ,x_ ,y_ ,u_ , v_ , sgn_ : -1] :=Block[{i,cm,sust,cev,lt}, ( 

cm=Complement [Union [Flatten [Outer [Plus , Apply [Dot , autono [1] , 1] , 

sgn Apply [Dot, 1,1] ] ] ] ,{0}] ; 

sust={a [10] ->x , a [20] ->y , a [30] ->u , a [40] ->v , 

a[13]->Adj [x] ,a[23]->Adj [y] ,a[33]->Adj [u] ,a[43]->Adj [v] , 

a [11] ->Transpose [x] ,a[31]->Transpose[u] , 

a [12] ->Conjugate [x] , a [32] ->Conjugate [u] } ; 

cmev=cm/ . sust ; 

For [lt={> ; i=l , i<=Length [cm] , i++ , 
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If [Chop [Det [cmev [ [i] ] ] ] ! =Conjugate [Chop [Det [cmev [ [i] ] ] ] ] , 
AppendTo[lt,cm[[i]]] ] ];lt)]; 

(* Definition of Newecs *) 

Newecs [oldecs_ , list_ ,x_ ,y_ ,u_, v_ , vars_ , complexvars_,ops ] : = 

Block [{ecs2 , ecs3 , i , j , It s , It , newec , end , what , detyn , trordet , vars2} , ( 
detyn= (det/. {ops}) /.Options [Newecs] ; 

If [detyn, trordet [args ] :=Imdet [args] , trordet [args ] : =Imtr [args] , 

trordet [args ] : =Imtr [args] ]; 

ecs2=oldecs ; 
ecs3={} ; 

vars2=Join [vars ,Re [complexvars] , Im [complexvars] ] ; 
For [lts=-Q ; i=l , i<=Length [list] , i++ , ( 
newec=trordet [list [ [i] ] ,x,y,u, v, complexvars] ; 
For [j=l;end=False;what=0, j<=Length[ecs2] && !end,j++, 
If [PolynomialQ [Factor [ecs2 [ [j] ] /newec] ,vars2] && ! (NumberQ [Factor 
[ecs2 [ [j] ] /newec] ] ) , what=j ; end=True] ] ; 
If [end, ecs2 [ [what] ] =newec ; AppendTo [Its, list [ [i] ] ] ] )] ; 
For [lt={> ; i=l , i<=Length [list] , i++ , ( 
newec=trordet [list [ [i] ] ,x,y,u,v, complexvars] ; 
For [j=l;end=False, j<=Length[ecs2] && !end,j++, 

If [PolynomialQ [Factor [newec/ ecs2 [ [ j ] ] ] , vars2] , end=True] ] ; 
If [! end, AppendTo [ecs3, newec] ; AppendTo [It , list [ [i] ] ] ] )] ; 
{ecs2,lts,ecs3,lt} )] ; 

Options [Newecs] ={det->False> 

(* Definition of Eliminatesols *) 

Eliminatesols [sols_] : =Block [{i , j , sols2 , inter} , ( 
sols2=Union [sols] ; 

For [i=l,i<=Length[sols2] ,i++,If [Head[sols2[[i]]]==And, , 
sols2 [ [i] ] =And [sols2 [ [i] ] , sols2 [ [i] ] ] , 
sols2 [ [i] ] =And [sols2 [ [i] ] , sols2 [ [i] ] ] ] ] ; 
i=l; 

While [i<Length[sols2] , ( 
j=i+l; 

While [j<=Length[sols2] , ( 
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inter=Intersection [sols2 [ [i] ] ,sols2[[j]] ] ; 
If [inter==Union [sols2 [ [i] ] ] , sols2=Delete [sols2 , j] , , 
If [inter==Union [sols2 [ [ j ] ] ] , 

sols2=Delete[sols2,i] ; j=i+l, , ] )] ; 

i++)]; 

For [i=l , i<=Length [sols2] , i++ , sols2 [ [i] ] =Union [sols2 [ [i] ] ] ] ; 
sols2 )] ; 

(* Definitions of Redl2 and Red21 *) 

Redl2 [m0_] :=Block[{lu,m,ml2,m21,m22,w2,w, delta, n>, ( 

lu={}; 

m=mO; 

ml2:={{m[[l,2]] ,m[[l,3]]}>; 
m21 : =Transpose [ml2] ; 

m22:={{m[[2,2]] ,m [ [2 ,3] ] } , {m [ [3,2] ] ,m[[3,3]]}}; 
w2=DiagonalizeS[m22] [[2]] ; 
w=IncDiml [w2] ; 
m=Chop [w . m . Transpose [w] ] ; 
If [Chop [m[[2,2]]-m[[3,3]]]==0, 
If [ml2=={{0,0}}, 

AppendTo [lu, 2] ; AppendTo [lu,3] , 
(* else *) 

If [m22=={{0,0},{0,0}}, 

w2=DiagonalizeH [m21 . Adj [m21] ] [ [2] ] ; 
w=IncDiml [w2] ; 
m=Chop [w . m . Transpose [w] ] , 
(* else *) 

If [m[[l,l]] !=0, 

delta— ArgCm[ [1,1]]] /2; 
w=DiagonalMatrix [{Exp [I delta] , 1 , 1}] ; 
m=w . m . Transpose [w] , 
(* else *) 
n= (Chop [ml2 . Transpose [ml2] ] ) [ [1 , 1] ] ; 
If [n==0, 

AppendTo [lu, 1] , 
(* else *) 

delta=-Arg [n] /2; 

w=DiagonalMatrix [{Exp [I delta] , 1 , 1>] ; 
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m=w . m . Transpose [w] 

] 

]; 

If [Im [ml2] =={{0,0}}, 

w2=DiagonalizeH[Re [m21] .Transpose [Re [m21] ] ] [[2]] ; 

w=IncDiml [w2] ; 

m=Chop [w . m . Transpose [w] ] , 
(* else *) 

w2=DiagonalizeH[Im[m21] .Transpose [Im[m21] ] ] [[2]] ; 

w=IncDiml [w2] ; 

m=Chop [w . m . Transpose [w] ] 

] 

] 

] (* end ml2=0 *) 
(* else 1st if *) ] ; 
m 

)]; 

Red21 [m_,me_] : =Block [{m2} , 
If [Chop [Abs [me [ [1 , 1] ] ] -Abs [me [ [2 , 2] ] ] ] ! =0 , 
m, 

(* else *) 

m2=Reverse [Transpose [Reverse [m] ] ] ; 

Reverse [Transpose [Reverse [Redl2 [m2] ] ] ] 
]]; 

(* Definition of LookCP *) 

LookCP [m_] : =Block [{m2 , dim , i , ecsdegO , ecsdegl , ecsdeg2 , allecs , sol} , ( 
m2=Chop [Arg [m] ] ; 
dim=Length [m] ; 

For [i=l,i<=dim,i++,If [m[[i,i]]==0,m2[[i,i]]=delta[i] ] ]; 
ecsdegO=Part [minors [m2] , 

Flatten [Position [Apply [Plus, Thread[vec [minors [m2] ,dim]] ,1] ,0]]] ; 
ecsdegl=Part [minors [m2] , 

Flatten [Position [Apply [Plus, Thread[vec [minors [m2] ,dim]] ,1] ,1]]] ; 
ecsdeg2=Part [minors [m2] , 

Flatten [Position [Apply [Plus, Thread[vec [minors [m2] ,dim]] ,1] ,2]]] ; 
allecs=Join [ecsdegl , ecsdeg2] ; 
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If [Complement [Chop [Sin [ecsdegO] ] , {0}] =={> , 

If [allecs=={} ,True , (* there aren't equations *) 

For [i=2 ; indeps={allecs [ [1] ] } , i<=dim && i<=Length [allecs] , i++ , 

If [linealindep [Append [indeps , allecs [ [i] ] ] , dim] , 

AppendTo [indeps , allecs [ [i] ] ] ] ]; 
sol=Flatten [Solve [indeps==Table [0, {Length [indeps] }] ] ]; 
Complement [Chop [Chop [Sin [allecs/ . sol] ] ] , {0}] =={>] , False] 
)] 



End[] 



EndPackage [] 
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